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INTRODUCTION 

A  constitutive  model  of  the  shock  compression  of 
snow  is  of  interest  for  its  direct  application  to  such  fields 
as  planetary  sciences,  cold  regions  and  military  engi¬ 
neering,  and  shock  isolation.  The  only  existing  data 
obtained  with  reliable  experimental  methods  are  the 
high-pressure  data  (3.8-35.4  GPa)  of  Bakanova  et  al. 
(1975).  Experiments  by  Napadensky  ( 1 964),  Wakahama 
and  Sato  (1977),  Sato  and  Brown  (1983)  and  Sato 
(1987)  did  not  meet  criteria  for  steady  plane- wave 
propagation  even  though  the  data  were  reduced  using 
that  assumption. 

We  have  conducted  a  test  program  using  embedded 
stress  gauges  to  obtain  stress-strain  relations  for  snow 
(Brown  et  al.  1988).  The  unsteady  nature  of  shock  wave 
propagation  in  the  snow,  the  specialized  sample  prepa¬ 
ration  methods,  and  the  large  impedance  mismatches 
between  the  snow  and  stress  gauges  resulted  in  complex 
stress  histories.  These  features  precluded  direct  appli¬ 
cation  of  the  Rankine-Hugoniot  jump  conditions  or  the 
Lagrangian  conservation  equations  for  mass,  momen¬ 
tum  and  energy  (Fowles  and  Williams  1970,  Seaman 
1976)  to  analyze  the  experimental  results.  Instead,  a 
detailed  analysis  of  one  of  our  experiments,  using  the 
PRONTO  2D  dynamic  finite-element  program  (Taylor 
and  Flanagan  1987),  was  conducted  to  determine  the 
origins  of  distinctive  features  from  measured  stress 
histories  and  to  determine  loading,  unloading  and  re¬ 
loading  paths  for  the  snow.  In  this  paper  we  discuss  our 
experimental  procedures  anddifficulties,  the  data  records 
for  a  typical  shot,  our  assumptions  for  determining  the 
initial  loading  curve,  and  unloading  and  reloading  curves 
for  the  snow,  as  well  as  our  results.  The  overall  shock 
response  of  snow  derived  from  the  entire  shot  test  series 
will  be  presented  separately. 


EXPERIMENTAL  METHODS 

A  shock  wave  experiment  was  conducted  using  a 
200-mm-diameter  gas  gun  to  obtain  uniaxial  strain 
loading  in  snow.  A  flat-faced  35.9-mm-thick,  polym¬ 
ethyl  methacrylate  (PMMA)  flyer  plate  was  accelerated 
with  compressed  nitrogen  gas  to  impact  a  snow  target. 

The  target  assembly  consisted  of  a  copper  cylinder 
with  sealed  aluminum  capping  end  plates  to  provide  a 
vacuum-tight  cannister  (Fig.  1).  A  vacuum-tight  target 
assembly  was  required  because  the  target  chamber  of 
the  gas  gun  was  evacuated  prior  to  firing.  The  snow, 
which  had  a  vapor  pressure  of  about  300  Pa,  would  have 
sublimated  if  it  had  not  been  isolated  in  a  sealed  cannis¬ 
ter.  A  spiral  of  copper  tubing  was  soldered  around  the 
copper  cylinder,  allowing  refrigeration  of  the  sample 
once  the  target  had  been  mounted  on  the  gas  gun.  The 
target  was  assembled  in  an  adjacent  coldroom  with  the 
target  axis  vertical.  A  buffer  plate,  consisting  of  a 
carbon  stress  gauge  epoxied  between  two  aluminum 
plates  (12.7-mm-thick  plate  epoxied  to  a  6.4-mm-thick 
plate),  formed  the  front  of  the  target. 

Snow  was  sieved  into  the  copper  cylinder  in  stages, 
with  stress  gauges  clamped  in  place  at  specific  distances 
from  the  aluminum/snow  interface  (Table  1 ).  The  stress 
gauges  were  placed  off-axis  so  that  the  shock  wave  had 
a  direct,  unobstructed  path  from  the  buffer  through  the 
snow  to  each  gauge.  Our  snow  had  approximately 
millimeter-sized  grains.  Consequently  our  uncertainty 
of  gauge  position  was  about  ±1  mm,  which  for  this 
experiment  corresponds  to  an  uncertainty  in  arrival  time 
of  about  ±5  ps. 

Thermocouples  were  placed  in  the  snow  and  were 
used  to  monitor  the  sample  temperature.  The  copper 
cylinder  was  completely  filled  with  snow,  and  the 
instrumentation  wires  were  run  along  the  inside  wall  of 


Figure  1.  Schematic  of  the  snow  target  as¬ 
sembly,  consisting  of  an  aluminum  plate  (A 
and  B),  a  carbon  gauge  (J),a  cooling  coil  (C ) 
soldered  to  a  copper  cylinder,  an  aluminum 
back  support  plate  (D),  gauge  and  thermo¬ 
couple  leads  exiting  the  rear  surface  (E),  and 
a  gauge  support  pedestal  mounted  to  the  inside 
wall  of  a  copper  cylinder  (F). 


the  copper  cylinder  and  through  vacuum-sealed  holes  in 
the  back  plate  (12.7-mm  aluminum).  Six  shorting  pins, 
each  with  a  different  length,  were  mounted  concentri¬ 
cally  on  the  outer  edge  of  the  aluminum.  These  shorting 
pins  were  used  to  determine  the  flyer  plate  impact 
velocity  and  to  trigger  data  acquisition. 

After  the  target  was  assembled,  the  sample  was  “cold 
soaked”  in  a  refrigerated  room  overnight  before  being 
mounted  on  the  end  of  the  gas-gun  barrel.  Snow  tem¬ 


perature  was  controlled  by  circulating  cold  nitrogen  gas 
through  the  copper  tubing  attached  to  the  outer  wall  of 
the  target. 

Stress-time  records  were  measured  using  50-ohm 
carbon-film  piezoresistive  gauges  (Krehl  1978).  Their 
high  sensitivity  allows  for  a  relatively  low  excitation 
power  and  less  Joule  heating  of  the  gauge  compared  to 
other  piezoresistive  gauges.  The  gauge  has  acalibration 
uncertainty  on  loading  and  unloading  of  between  5  and 


Table  1.  Stress  gauge  positions  relative  to  the  snow/aluminum  interface. 


Gauge 

plane 

Gauge 

number 

Relative 

position 

(mm) 

Gauge 

thickness 

(mm) 

Average  relative 
position 
(mm) 

0 

1 

-6.45  ±0.01 

0.1778* 

-6.45 

1 

2 

13.2  ±  1 

0.8128 

14.0 

3 

14.8+  1 

0.6604 

2 

4 

27.2  ±  1 

0.6858 

26.0 

5 

24.8  ±  1 

0.5080 

3 

6 

37.2  ±  1 

0.5824 

37.2 

‘Includes  epoxy  and  mylar  tape  (the  unclad  thickness  was  0.0762  mm). 
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10%  of  the  stress  at  stresses  below  about  2.0  GPa 
(Gourdin  and  Weinland  1986,  King  and  Janie  1987). 
Hysteresis  effects  range  from  0  to  9%  of  the  stress  on 
partial  or  complete  unloading  (King  and  Janie  1987). 
Our  experience  in  using  the  gauge  in  soils  and  alumi¬ 
num  have  not  shown  evidence  of  measurable  hysteresis. 
The  active  element  for  the  carbon  gauge  (0.75  x  1 .25  cm) 
forms  a  single,  continuous,  wide  strip  rather  than  a  grid, 
as  is  the  case  for  the  more  commonly  used  manganin 
and  ytterbium  gauges;  it  is  thus  less  susceptible  to 
puncturing  by  individual  snow  grains.  The  gauges  were 
encapsulated  between  0.025-;  im-thick  layers  ofkapton. 
Recording  life  was  extended  by  using  a  0.19-  to  0.34- 
mm  layer  of  mica  as  armor. 

Pulsed  Wheatstone  bridge  power  supplies  were  used 
to  provide  70  V  of  gauge  excitation  for  200  ps.  Bridge 


output  signals  were  recorded  digitally  using  CAMAC- 
based  waveform  digitizers  with  a  band  width  of  1  MHz 
and  a  sampling  rate  of  2  MHz. 


EXPERIMENTAL  RESULTS 
AND  DIFFICULTIES 

Figure  2  shows  data  from  the  six  stress  gauges  in  the 
target  assembly  (Fig.  1 ):  one  in  the  aluminum  buffer  and 
the  remainder  in  the  snow  (compressive  stresses  are 
treated  as  positive).  For  this  shot  the  snow  had  an  initial 
density  of  400  kg  m-3  and  a  temperature  of  -8°C.  The 
flyer  plate  impact  velocity  was  150.7  ms-1  .resulting  in 
about  a  0.4-GPa  peak  impact  stress  with  four  reverbera¬ 
tions  in  the  target  buffer  (Fig.  2,  gauge  1  at  gauge  plane 


Time  (ps) 


Figure  2.  Stress  records  from  gauge  1  in  the  aluminum  and  gauges 
2-6  embedded  in  snow. 
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0).  Shock  propagation  velocities  in  the  snow  were  242 
±  15  m  s_l  at  gauge  plane  1  (gauges  2  and  3),  230  ±  15 
m  s_1  at  gauge  plane  2  (gauges  4  and  5)  and  205  ±  15  m 
s-1  at  gauge  plane  3  (gauge  6).  Gauge  records  from  the 
first  gauge  plane  show  an  initial  spike  at  shock  arrival, 
a  gradual  decrease  to  a  minimum  value  of  about  7-8 
MPa  followed  by  a  shallow  rise  to  about  8-10  MPa.  A 
slight  dip  in  stress  level  then  occurs,  with  a  spike  at  1 35 
jj.s  before  a  decrease  in  stress  below  7  MPa.  The  stress 
records  at  the  second  and  third  gauge  planes  show  a 
significantly  longer  initial  rise  time  and  smaller  spikes. 
Rankine-Hugoniot  theory  predicts  a  stress  level  of 
about  12  MPa  for  a  steady  wave  traveling  at  240  m  s-> 
in  the  snow  if  the  snow  particle  velocity  was  equal  to  the 
maximum  possible  buffer  particle  velocity  after  impact 
with  the  PMMA  flyer  plate  (about  120  m  s-1). 

The  generally  good  agreement,  within  experimental 
uncertainty,  between  the  stress  records  in  the  same 
plane  for  gauges  2  and  3,  and  gauges  4  and  5 ,  and  the  fact 
that  pseudo-steady  stresses  for  gauges  2  and  3  are  on  the 
same  order  as  estimated  using  Rankine-Hugoniot  theory, 
suggest  that  the  measured  records  are  valid. 

The  stress  records  from  all  of  the  gauges  in  the  snow 
are  both  complex  and  unsteady.  Two  factors — target 
design  and  impedance  mismatching  between  the  snow 
and  stress  gauges — contribute  to  the  complexity.  Our 
target  design  was  constrained  by  the  need  fora  vacuum- 
tight  container  that  would  still  produce  a  one-dimensional 
strain  shock  when  impacted.  We  achieved  a  flat  impact 
surface  by  using  a  thick  buffer  of  sufficient  strength  to 
prevent  any  significant  deflections  due  to  the  pressure 
difference  between  the  interior  of  the  target  and  the 
target  chamber.  Upon  impact  from  the  flyer  plate,  the 
buffer  imparted  stresses  to  the  snow,  which  gradually 
decreased  as  the  momentum  was  transferred  from  the 
buffer  to  the  target  by  multiple  reverberations  in  the 
buffer.  The  rate  of  decrease  in  applied  stress  was  con¬ 
trolled  by  the  impedance  mismatch  between  the  buffer 
and  the  snow.  A  large  impedance  mismatch  between  the 
stress  gauges  and  the  snow  produced  large-amplitude 
reflected  pulses,  furthercomplicating  the  signal  (Brown 
et  al.  1988,  Gaffney  1989).  Shock  wave  stress  records 
were  sufficiently  long  that  waves  generated  at  the  edge 
of  the  target  propagated  to  the  stress  gauges  prior  to  the 
end  of  the  experiment. 


MODEL  SIMULATION  OF  THE 
EXPERIMENTAL  DESIGN 

We  used  the  PRONTO  2D  transient  solid  dynamics 
finite-element  program  as  a  tool  to  understand  the 
complex  wave  records  and  to  construct  constitutive 
models  of  shock  propagation  in  snow.  We  initially 


constructed  a  simple  model  incorporating  only  the 
flyer,  buffer  and  snow  in  a  one-dimensional  ( 1 D)  geom¬ 
etry.  Ultimately  it  was  necessary  to  include  the  stress 
gauges  and  the  failure  strength  of  the  epoxy  bond 
between  the  two  plates  that  made  up  the  buffer  in  the  1 D 
model.  Figure  3  shows  the  final  model  geometry  used  in 
our  1 D  analysis.  In  the  model  a  PMMA  flyer  impacted 
an  aluminum  buffer  that  had  a  no-friction  contact  with 
the  snow.  The  aluminum  buffer,  in  the  model,  consisted 
of  a  mylar  layer  sandwiched  between  two  plates  of 
aluminum  tooling  plate,  which  simulated  the  epoxy, 
mylar  tape  and  stress  gauge  between  the  two  buffer 
plates  used  in  the  test  experiment.  Mica-clad  gauges 
were  simulated  by  layers  of  mica  at  the  average  position 
of  gauges  2  and  3  for  the  first  gauge  plane  and  gauges  3 
and  4  for  the  second  gauge  plane.  The  third  gauge  plane 
was  placed  at  the  same  position  in  the  finite-element 
mesh  as  the  gauge  6  position  in  the  actual  test.  Gauge 
dimensions  used  in  the  model  were  taken  from  gauges 
2, 4  and  6  (Table  1). 

The  PMMA  and  aluminum  were  modeled  as  an 
elastic-plastic-hydrodynamic  material.  This  model 
combines  a  Von  Mises  yield  condition  including  strain 
hardening  to  describe  the  deviatoric  response.  Volu¬ 
metric  response  is  represented  by  a  Mie-Gruneisen 
equation  of  state  (EOS)  of  the  form 

p  =  *o  n  ( i  +  *1  ti  +  k2  if )  ( i  -  Eo. )  +  £v 

2  P 

(1) 

where 

£v  =  energy  per  unit  volume 

p0  =  initial  density 
p  =  density  at  pressure  P 
I)  =  (l-po/p) 

and  the  remaining  parameters  are  material-dependent 
properties  and  are  given  in  Table  2.  Mylar  was  modeled 
as  a  pure  hydrodynamic  material  using  the  EOS  param¬ 
eters  tabulated.  Mica  was  treated  as  an  elastic  material. 
Although  both  deviatoric  behavior  and  volumetric  be¬ 
havior  of  the  materials  were  included  in  the  model, 
changes  in  deviatoric  parameters  had  a  negligible  effect 
on  the  calculated  results. 

Modulation  of  the  uniaxial  strain  conditions  in  our 
impact  experiment  by  release  waves,  generated  at  the 
lateral  edge  of  the  target,  was  examined  using  an 
axisymmetric  two-dimensional  (2D)  model  (Appendix 
A).  Material  parameters  and  the  constitutive  model  for 
snow  were  the  same  for  both  the  2D  and  1 D  models. 

In  constructing  our  model  we  used  the  measured 
stresses  as  a  guide  to  establish  the  validity  of  calculated 
results.  The  measured  stress-time  record  for  gauge  1 
(Fig.  2)  was  used  to  de^imine  if  the  calculated  stresses 
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Figure  3.  Schematic  of  the  experimental  configuration  used  to  generate  the  one-dimensional  finite- 
element  mesh  for  PRONTO  2D.  Gauge  planes  0,1,2  and  3  are  designated  hy  GPO,  GP1 .  GP2  and 
GP3  in  the  figure.  The  mesh  was  one  element  wide,  with  the  number  of  lateral  elements  for  each 
component  shown  in  the  figure. 


from  the  model  produced  by  the  impact  of  the  PMMA 
flyer  and  aluminum  buffer  were  realistic.  Figure  4 
shows  a  comparison  of  the  model-calculated  stresses  in 
the  buffer,  at  the  location  of  gauge  1 ,  and  the  measured 
stresses  determined  from  gauge  1 .  Calculated  stresses 
in  the  buffer  were  not  significantly  influenced  by  the 
snow  behavior  for  any  reasonable  choice  of  snow  pres¬ 
sure  vs  volumetric  strain  (P~\r)  curve.  The  calculated 


and  measured  stresses  agree  very  well  for  the  first  peak 
and  fall.  Thereafter  there  is  fair  agreement  in  the  timing 
of  the  respective  maxima  and  minima.  The  four  rever¬ 
berations  calculated  by  the  model  are  caused  by  re¬ 
flected  waves  from  the  aluminum/snow  interface  and 
the  PMMA/aluminum  interface.  The  reverberations 
end  when  the  unloading  wave  from  the  back  of  the 
PMMA  flyer  arrives  at  the  mylar  layer,  causing 


Table  2  .  Constitutive  model  parameters  used  to  develop  the  PRONTO  2D  model  calculation. 


Young's  Yield  Hardening  Pressure 


Material 

Model  type 

Density 
(kgm ->) 

modulus 

(GPa) 

Poisson  s 
ratio 

strength 

(MPa) 

modulus 

(MPa) 

Bela 

cutoff 

(MPa) 

k» 

(GPa) 

k, 

Go 

PMMA* 

Elastic/plastic- 

1184 

6.3 

0.4 

375 

500 

0.5 

-100 

5.83 

4.42 

13.04 

0.7 

hydrodynamic 

Aluminumt 

Elastic/plastic- 

2828 

72.5 

0.33 

270 

1150 

0.5 

-1000 

81.86 

2.81 

7.25 

2 

hydrodynamic 

Mylar** 

Hydrodynamic 

1390 

6.73 

1.196 

2.29 

1 

Micatt 

Elastic 

2844 

69.6 

0.26 

•PMMA  EOS  data  (Rice  1980)  were  reformulated  to  conform  to  the  PRONTO  2D  representation  of  the  Mie-Gruneisen  EOS.  and  the  hardening 
modulus  was  estimated.  Beta  is  a  parameter  that  determines  the  amount  of  isotropic  and  kinematic  hardening  and  was  estimated. 
tAluminum  EOS  from  Rice  (1980).  The  hardening  modulus  was  determined  from  data  given  by  Herrmann  and  Lawrence  (1978).  Beta  was 
estimated. 

••Mylar  EOS  from  Louie  et  al.  (1970). 

ttCalculated  from  isotropic  elastic  moduli  (Vaughan  and  Guggenheim  1986). 
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Figure  4.  Comparison  of  the  calculated  stresses  for  gauge  plane  0  with  the  measured 
stress  record  for  gauge  l,  located  in  the  aluminum  buffer. 


debonding.  Calculations  have  shown  that  tensile  stresses 
in  the  buffer  would  have  exceeded  120  MPa  had  the 
band  between  the  mylar  and  the  aluminum  not  been 
allowed  to  fail.  The  reduction  in  calculated  peak  stress 
for  each  succeeding  reverberation  is  controlled  by  the 
impedance  mismatches  at  the  aluminum/PMMA  flyer 
interface  and  the  aluminum/snow  interface.  The  most 
probable  explanation  for  the  difference  between  mea¬ 
sured  and  calculated  stresses  for  the  later  reverberations 
is  that  gauge  1  and/or  the  epoxy/mylar  bond  was  dam¬ 
aged  after  being  subjected  to  a  tensile  stress  of  about  25 
MPa  at  4.8  ps  (Fig.  4).  Also  the  carbon  gauges  are  not 
designed  to  measure  tensile  stresses.  Epoxy  resins,  such 
as  we  used  in  manufacturing  the  buffer,  have  tensile 
strengths  that  range  from  14  to  90  MPa,  depending  on 
their  composition  (Baumeister  and  Marks  1967).  After 
the  first  tensile  pulse  the  gauge  was  still  able  to  sense 
compressive  stresses  but  not  tensile  stresses.  The  gauge 
failed  completely  upon  the  arrival  of  the  unloading 
wave  from  the  flyer  (Fig.  4).  This  unloading  wave  also 
caused  the  two  pans  of  the  buffer  plate  to  separate.  The 
low-amplitude  stresses  (less  than  10  MPa)  in  the  calcu¬ 
lated  record  after  the  arrival  of  the  unloading  wave, 
shown  in  Figure  4,  are  caused  by  reverberations  in 
buffer  disk  B  (in  contact  with  the  snow,  Fig.  1 ).  These 
reverberations  cause  the  gradual  decrease  in  stresses  in 
the  snow  after  shock  front  passage. 

Measured  stress-time  records  from  gauges  embed¬ 
ded  in  the  snow  were  used  to  aid  in  our  construction  of 
a  constitutive  relationship  for  snow  compaction.  This 
was  accomplished  by  assuming  that  snow  strain-hard- 


ens  during  all  stages  of  its  compaction  and  that  the 
form  of  the  dynamic  P-V  curve  is  similar  to  that  found 
from  quasi-static  uniaxial  strain  tests  (Abele  and  Gow 
1975).  That  is,  the  slope  of  the  stress-strain  curve  is 
either  increasing  or  constant  with  increasing  volumetric 
strain. 

We  used  the  Soils  and  Crushable  Foams  material 
model  in  PRONTO  2D  to  describe  snow  compaction. 
This  phenomenological  model  was  constructed  from  a 
version  of  the  model  developed  by  Krieg  (1978)  and 
described  by  Swenson  and  Taylor  ( 1 983).  It  consists  of 
a  yield  function  describing  deviatoric  deformations 
uncoupled  from  the  P-V  curve  used  to  define  volumet¬ 
ric  deformations.  The  yield  surface  oy  is  a  surface  of 
revolution  about  the  hydrostat  specified  as  a  quadratic 
in  pressure 

Oy  =  <^)  +  a|  P  +  ai  P*  .  (2) 

We  used  ao  =  10  kPa,  a\  =0.1  and  02  =  0.0  in  our  cal¬ 
culations.  The  choice  of  parameters  in  the  yield  function 
does  not  greatly  affect  the  results  since  the  deformations 
are  dominated  by  volumetric  compaction.  Allowing  the 
yield  function  to  depend  slightly  on  pressure  does, 
however,  greatly  increase  the  stability  of  the  calcula¬ 
tion. 

The  user-defined  function  /p(£(P*x  ).  where  the  vol¬ 
umetric  strain  is  ev  =  -  In  ( p^p),  is  used  to  define  the 
P-V  curve  (Taylor  and  Flanagan  1987).  This  function 
incorporates  both  a  loading  curve  and  a  linear  unload¬ 
ing-reloading  curve.  The  slope  of  the  unloading  curve 
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is  required  to  be  greater  than  the  loading  P-V curve.  The 
function 

<t)p=/>-/p(erx)  (3) 

where  6™“  is  the  maximum  ev  experienced  by  each 
element,  defines  the  motion  of  the  end  cap  of  the 
deviatoric  yield  function  along  the  hydrostat.  Tensile 
failure  is  assumed  to  occur  if  the  pressure  is  less  than  a 
user-defined  fracture  pressure  Pfr.  Tensile  failure  strain 
£fr  is  initially  negative  and  is  set  to  Pfr/Ko,  where  Kq  is 
the  initial  bulk  modulus  of  the  snow;  £fr  increases  with 
E^**  and  is  determined  by  unloading  from  the  /p(£™ax  ) 
curve  along  a  linear  unloading-reloading  curve  to  PfT. 
The  strain  at  Pfr  is  the  new  value  of  £ff. 

The  results  of  preliminary  calculations  using  the 
Soils  and  Crushable  Foams  model  replicated  many 
features  of  our  measured  data,  including  arrival  times 
and  peak  stresses  that  are  in  good  agreement  with 
measurements  (Johnson  et  al.  1990).  However,  the 
post-peak  stress  histories  could  not  be  reproduced  be¬ 
cause  of  the  restrictions  that  the  unloading-reloading 
curve  be  linear,  have  the  same  value  as  the  initial  bulk 
modulus  irrespective  of  the  total  volumetric  strain,  and 
have  a  slope  greater  than  the  maximum  slope  of  the  P- 
V  curve.  We  modified  the  Soils  and  Crushable  Foams 
model  to  accommodate  more  general  unloading  and 
reloading  behavior  for  snow.  On  initial  loading,  the 
pressure  in  the  snow  is  a  function  of  E™1*  ,  as  was  the 
case  in  the  original  model.  Multiple  linear  or  nonlinear 
unloading  paths  may  be  defined,  with  each  unloading 
path  being  defined  over  a  specified  range  of  £v.  One  of 
two  possible  reloading  schemes  may  be  used:  reloading 
along  the  same  path  as  unloading  or  reloading  along  a 
chord  from  the  last  P-V  state  to  the  /p(£(J’ax  ),  e™3* 
state  (Appendix  B). 

CONSTITUTIVE  RELATIONSHIP 
AND  DISCUSSION 

Development  of  a  constitutive  relationship  proceeded 
in  two  stages  by  first  finding  a  P-V  curve  that  produced 
a  good  agreement  between  calculated  and  measured 
shock  arrival  times.  Next,  le  stress  history  after  shock 
arrival  was  used  to  constrain  the  unloading  and  reload¬ 
ing  curves.  The  P-V  curve  and  the  unloading  and  re¬ 
loading  curves  found  by  comparing  calculated  and 
measured  results  are  not  unique.  They  can  vary  because 
of  experimental  uncertainties,  the  complex  relationship 
of  one  segment  of  the  P-V  curve  to  another,  and  the 
subjective  determination  of  what  constitutes  a  good 
model  simulation  of  the  measured  results. 

The  range  of  possible  P-V  curve  shapes  was  con¬ 


strained  in  this  study  by  several  factors.  We  assumed 
that  the  slope  of  the  P-V  curve  either  increases  or  re¬ 
mains  constant  with  increasing  volumetric  strain  (which 
is  consistent  with  observations  from  quasi-static  tests 
on  snow  assuming  no  phase  changes  occur).  We  have 
included  the  further  constraint  that  the  final  density  for 
the  snow  must  be  less  than  that  of  our  recovered  post- 
experiment  samples  (860  kg  nrr3).  In  addition,  our 
calculations  have  shown  that  the  P-V  curve  can  be 
broken  into  two  regions,  which  affect  the  shock  wave 
differently.  The  magnitude  of  the  P-V' curve  in  the  mid¬ 
stress  range  (2.5  MPa  <  P  <  16  MPa)  affects  the  shock 
velocity  and  produces  moderate  changes  in  the  magni¬ 
tude  of  post-peak  stresses.  The  magnitude  of  the  P-V 
curve  in  the  high-stress  range  (16  MPa  <  P  <  40  MPa) 
does  not  affect  shock  velocity  but  is  important  in  re¬ 
producing  the  shock  arrival  stress  spikes  at  the  stress 
gauges  and  does  affect  the  post-peak  stress  amplitudes. 
The  P-V  curve  definition  in  the  low-stress  range  could 
not  be  determined  directly  from  our  tests.  The  curves 
were  estimated  from  the  quasi-static  test  results  of 
Abele  and  Gow  (1975).  Unloading  and  reloading  curves 
were  constrained  by  the  shock  wave  arrival  times  and 
stress  magnitudes  of  waves  reflected  from  gauges  in  the 
snow  and  from  the  alunrnum/snow  interface.  The  un¬ 
loading  and  reloading  curves  were  further  constrained 
by  comparing  the  attenuation  of  the  calculated  peak 
stresses  and  the  magnitude  and  frequency  of  post-peak 
stress  oscillations  with  those  measured. 

Initial  calculations  used  a  P-V  curve  from  a  quasi¬ 
static  compaction  test  by  Abele  and  Gow  (1975)  (Table 
3).  These  produced  arrival  times  that  were  much  later 
than  measured  for  all  gauges  in  the  snow.  Succeeding 
calculations  were  made  after  increasing  the  model  P-V 
curve  slope  in  a  systematic  manner  (that  is,  increasing 
the  stress  needed  to  cause  a  given  volumetric  strain). 
The  final  P-V  curve  and  the  unloading  curves  used  to 
calculate  the  results  presented  in  this  paper  are  shown  in 
Figure  5  and  listed  in  Table  4.  The  calculated  peak 
stresses  for  gauges  2, 3, 4  and  5  and  the  post-peak  stress 
history  for  gauge  1  agree  with  measured  values  within 
the  limits  of  experimental  error.  The  calculated  stresses 
for  gauge  6  are  larger  than  for  the  measured  record  but 
have  the  same  form.  Differences  between  the  calculated 
stress  histories  for  gauges  4, 5  and  6  and  the  measured 
values  may  result  from  superposition  of  release  waves 
that  originated  at  the  edge  of  the  aluminum,  degrading 
the  uniaxial  strain  conditions. 

Comparison  of  the  ID  calculated  and  measured 
stresses  for  gauges  2  and  3  are  shown  in  Figure  6.  This 
provides  a  means  of  identifying  important  processes 
affecting  the  shock  propagation.  Calculations  show  that 
impedance  mismatch  between  the  gauges  and  the  snow 
produced  the  stress  spike  at  shock  arrival.  The  measured 
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Table  3.  P-V  curve  parameters  for  Abele  and  Gow’s  measurements  (1975). 


fp(£  vac)  Density  2\s  K0  Py, 

(MPa)  ( kg  in 3'  ty,  (MPa)  (MPa)  (MPa) 


0.0 

0.0 

0.109 

0.14 

461 

0.167 

0.28 

570 

0.379 

0.7 

581 

0.399 

1.4 

649 

0.509 

2.1 

680 

0.556 

2.88 

719 

0.612 

•The  shear  modulus  p  was  estimated  from  acoustic  data  on  snow. 

tThe  initial  bulk  modulus  K0  was  chosen  to  be  about  the  same  as  the  bulk  modulus  of  ice. 


stress  after  shock  arrival  decreases  to  that  of  the  ambient 
snow  stress.  Measured  stresses  then  increase  at  95  (is  and 
remain  elevated  until  about  120  (is.  Calculations  indi¬ 
cate  that  the  stress  rise  from  95  to  1 20  |is  is  due  to  a  shock 
pulse  that  had  been  reflected  from  gauge  plane  1  propa¬ 
gated  back  to  the  aluminum/snow  interface  and  reflected 
back  to  gauge  plane  1 .  The  doubly  reflected  wave  shows 
the  effects  of  dispersion  and  attenuation  during  its 
travel.  At  1 20  ps  the  measured  stress  begins  to  decrease 
because  of  unloading  from  the  buffer.  The  measured 


stress  decrease  is  interrupted  by  the  arrival  of  a  shock 
pulse  that  calculations  indicate  originated  from  a  re¬ 
flected  wave  at  the  second  gauge  plane  (the  gauge 
located  at  0.0514  m  in  Fig.  3,  gauge  plane  2).  This  pro¬ 
duces  a  sharp  stress  rise  at  135  (is,  then  a  stress  de¬ 
crease  to  the  stress  magnitude  of  the  buffer  unloading 
wave. 

Calculations  using  the  2D  model  (Appendix  A) 
show  that  gauge  plane  1  stresses  are  not  significantly 
affected  by  release  waves  originating  at  the  edge  of  the 


Figure  5.  Loading  pressure-volumetric  strain  curve  and  unloading¬ 
reloading  curves  used  in  PRONTO  2D  to  model  snow  compression  and 
Abele  and  Gow’s  quasi-static  pressure-volumetric  strain  results. 
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Table  4.  Snow  P-V  curve  parameters,  moduli  data  and  unloading  curve  parameters  used  to  define  snow  behavior  in 
PRONTO  2D  model  calculations. 


5>|elJaXJ  f(i](Ev)  fic{6v) 


Density 
(kg  nr3) 

P* 

(MPa) 

Pm, 

(MPa) 

K/m 

(MPa) 

Pirn 

(MPa) 

(MPa) 

Pirn 

(MPa) 

Kjm 

(MPa) 

2\l 

(MPa) 

Ko 

(MPa) 

Pfr 

(MPa) 

400 

0 

0.0 

10 

55 

0.681 

15 

217 

0.718 

15 

800 

0.764 

710 

4.2 

-0.01 

550 

0.35 

0.318 

15 

80 

20 

265 

30 

835 

635 

5.0 

0.462 

40 

2700 

740  10.0  0.615 

790  15.0  0.681 

820  20.0  0.718 

840  25.0  0.742 

848  30.0  0.751 

859  40.0  0.764 


aluminum  until  about  30  (is  after  the  peak  stress  arrival. 
Stresses  are  reduced  for  the  following  25  ps,  then  in¬ 
crease  to  slightly  larger  magnitudes  than  expected  for 
ID  wave  propagation.  This  is  consistent  with  the 
measured  stresses  for  gauges  2  and  3,  which  show 
pronounced  stress  reductions  about  20  ps  after  peak  stress 
arrival.  Stresses  increase  again  20-30 ps  later,  just  prior 
to  the  arrival  of  the  reflected  pulse  from  gauge  plane  2. 
We  cannot  unambiguously  interpret  the  measured  stress 
rise  for  gauges  2  and  3  at  about  1 10  ps  since  ID  model 
calculations  indicate  that  reflected  waves  from  the 
aluminum/snow  interface  arrive  at  between  95  and  1 20 
ps,  overlapping  the  arrival  of  the  2D  release  waves. 


Figure  7  compares  the  calculated  stress  history  for 
gauge  plane  2  with  measured  stresses  from  gauges  4  and 
5.  The  peak  calculated  stress  is  about  the  same  as  was 
measured.  The  2D  analysis  indicates  that  release  waves 
begin  affecting  the  shock  wave  upon  arrival  at  gauge 
plane  2.  Differences  between  the  post-peak  calculated 
and  measured  stresses  are  attributed  to  release  wave 
superposition  on  the  ID  stresses. 

Calculated  stresses  for  gauge  plane  3  have  the  same 
form  but  larger  amplitude  than  the  stresses  measured  by 
gauge  6  (Fig.  8).  Analysis  using  the  2D  model  indicates 
that  release  waves  propagating  from  the  aluminum 
buffer  reduce  stress  magnitudes  at  gauge  plane  6  but  do 


Figure  6.  Comparison  of  the  calculated  stresses  for  gauge  plane  1  with  measured  stress 
records  for  gauges  2  and  3. 
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Figure  7.  Comparison  of  the  calculated  stresses  for  gauge  plane  2  with  measured 
stress  records  for  gauges  4  and  5. 


not  affect  the  wave  form  significantly.  This  suggests 
that  magnitude  differences  between  measurements  and 
calculations  arc  caused  by  the  2D  edge  effects. 


CONCLUSION 

The  complex  unsteady  stress  histories  from  a  uniaxial 
strain  impact  on  snow  have  been  simulated  using  the 
PRONTO  2D  transient  solid  dynamics  finite-element 
program. 

The  P-V  curve  estimated  using  PRONTO  2D  is  not 


unique  and  can  vary  because  of  experimental  uncertain¬ 
ties,  the  interdependent  influences  of  different  regions 
of  the  curve,  and  the  subjective  determination  of  what 
constitutes  a  good  model  simulation.  The  range  of 
possible  P-V  curve  shapes  is  constrained  by  separate 
physical  properties  in  two  of  the  major  stress  regions. 
The  mid-stress  range  of  the  curve  is  determined  by 
shock  velocity  and  pseudo-steady  stress  levels.  The 
high-stress  region  of  the  curve  is  constrained  by  the 
magnitudes  of  stress  spikes  at  the  gauges.  The  low- 
stress  region  is  estimated  from  the  quasi-static  results 
of  Abele  and  Gow  (1975).The  shape  of  the  P-V  curve 


Figure  8.  Comparison  of  the  calculated  stresses for  gauge  plane  3  with  measured  stress 
records  for  gauge  6. 
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is  further  constrained  by  the  assumptions  that  slope 
either  increases  or  remains  constant  with  increasing 
volumetric  strain  and  that  the  final  density  for  the  snow 
must  be  less  than  860  kg  m-3.  Calculations  have  shown 
that  slight  variations  in  P-V  curve  shape  do  not  sig¬ 
nificantly  change  the  calculated  stress  histories.  Abele 
and  Gow’s  (1975)  quasi-static  P-V  curve  lies  well 
below  the  range  consistent  with  the  shock  data,  imply¬ 
ing  a  strong  strain-rate  dependence  for  shock  wave 
deformation  in  snow.  Post-peak  stresses  are  strongly 
controlled  by  the  unloading  and  reloading  behavior  of 
the  snow.  The  unloading  curve  for  snow  is  a  nonlinear 
function  of  volumetric  strain.  Reloading  is  along  a 
chord  from  the  last  P-V  state  to  the  maximum  pressure 
[/p(£viax  )]  and  maximum  volumetric  strain  (£™ax) 
state  (Fig.  Bl)  producing  a  hysteretic  unloading-re¬ 
loading  cycle  that  is  consistent  with  quasi-static  unload¬ 
ing  and  reloading  in  other  porous  media  (Johnson  and 
Green  1976,  Larson  and  Anderson  1979). 
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APPENDIX  A:  Propagation  of  Stress  Waves  from  the  Lateral  Edges  of  the  Aluminum 
Buffer  and  Snow 

OurPMMA  flyer  and  target  assembly  was  designed  to  have  a  relatively  large  diameter 
(about  200  mm)  so  that  release  waves  generated  at  the  lateral  edge  of  the  target  would  not 
arrive  at  the  stress  gauge  locations  until  after  the  experiment  ended.  Stress  histories  recorded 
by  the  gauges  in  our  experiment  were  sufficiently  long  that  distortion  of  the  uniaxial  strain 
conditions  by  release  waves  before  the  experiment  ended  was  likely.  We  examined  the 
propagation  of  release  waves  from  the  edges  of  the  target  assembly  into  the  center  axis  of  the 
snow  sample  by  constructing  an  axisymmetric  2D  model  of  our  experiment.  The  2D  model 
was  about  260  mm  in  diameter,  which  corresponded  to  the  diameter  of  the  buffer  plate  (a 
radius  of  about  1 30  mm  to  the  axis  of  symmetry),  and  had  the  same  configuration  on  axis  (with 
the  exception  that  no  gauges  were  embedded  into  the  snow),  material  properties  and 
dimensions  as  the  ID  model  (Fig.  Al). 

The  main  objective  of  conducting  the  2D  modeling  was  to  develop  a  qualitative 
understanding  of  how  release  waves  could  distort  the  stresses  in  our  shock  impact  experi¬ 
ment.  In  our  discussion  of  the  2D  modeling  results  we  discussed  pressure  rather  than  the 
stresses  parallel  and  perpendicular  to  the  center  axis  of  the  sample  since  those  stresses  are 
about  equal  to  each  other  and  to  the  pressure. 

The  pressure  distribution  from  the  edge  of  the  target  to  its  center  axis  as  a  function  of  time 
for  gauge  plane  0  located  in  the  aluminum  buffer  is  shown  in  Figure  A2a.  Pressure 
distributions  for  the  first  snow  element  adjacent  to  the  aluminum  buffer  and  at  the  location 
of  the  gauge  plane  1  in  the  snow  are  shown  in  Figures  A2b  and  A2c.  Figure  A2a  shows  the 
four  compressive  pulses  (also  seen  in  the  ID  model)  associated  with  the  initial  impact  of  the 
PMMA  flyer  with  the  aluminum  buffer  and  subsequent  reflections  from  the  aluminum/snow 
interface  and  PMMA/aluminum  interface.  The  influence  of  the  free  edge  of  the  aluminum 
buffer  can  be  clearly  seen  as  the  stress  pulses  are  reduced  near  the  free  edge  and  release  waves 
move  toward  the  center  axis  with  time.  The  diagonally  oriented  trough  in  Figure  A2a  and  the 
pressure  ridge  (A)  are  the  result  of  the  release  wave  propagating  towards  the  center  axis  from 
the  free  surface  with  a  propagation  velocity  of  4.0  ±  0.5  km  s_1 .  The  release  wave  reaches  the 


Figure  Al .  Schematic  of  the  two-dimensional  axisymmetric finite-element  mesh for  PRO  NT O 
2D.  Gauge  plane  0  is  designated  by  GP0  in  the  figure. 
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center  at  approximately  46  (is.  The  pressure  enhancement  at  the  center  axis  (A)  results  from 
superposition  of  the  corresponding  release  wave  from  the  opposite  side  of  the  aluminum 
buffer.  The  several  release  waves  seen  in  Figure  A2a  (A',  A,  B,  C)  can  be  traced  back  to  their 
origin  at  the  free  edge  of  the  aluminum  buffer  and  are  associated  with  the  four  primary 
pressure  pulses  caused  by  the  flyer’s  impact  with  the  aluminum  buffer. 

Figure  A2b  shows  that  the  release  waves  generated  in  the  aluminum  buffer  are  transmitted 
across  the  aluminum/snow  interface.  The  pressure  profile  of  the  release  waves  in  the  snow 
is  similar  to  that  of  the  release  waves  in  the  aluminum  buffer  (see  common  features  A,  B,  C 
and  D  in  Fig.  A2a  and  A2b).  Release  wave  pressures  are  greatly  attenuated  when  transmitted 
into  the  snow.  The  release  wave  generated  at  the  lateral  surface  of  the  snow  at  feature  E  (Fig. 
A2b)  is  propagating  inward  at  about  300  m  s_1.  At  this  velocity  the  release  wave  does  not 
reach  the  center  axis  until  about  450  |is,  well  after  the  experiment  has  ended. 

Figure  A2c  shows  that  the  gross  features  from  Figure  A2b  are  still  evident,  such  as  the 
trough  D;  however,  they  are  smoothed  and  attenuated  considerably.  The  release  wave 
generated  at  the  snow’s  lateral  edge  F  can  also  be  seen  propagating  toward  the  center  axis. 

Our  analysis  indicates  that  release  waves  originating  in  the  aluminum  buffer  arrive  at  the 
stress  gauge  planes  embedded  in  the  snow  before  the  experiment  ends.  At  gauge  plane  1  the 
release  waves  arrive  at  about  96  |0.s  and  modulate  the  stress  signal  by  about  1  MPa  (the  shock 
wave  arrives  at  gauge  plane  1  at  65  (is).  The  release  waves  begin  modulating  the  shock  wave 
front  before  it  arrives  at  either  gauge  plane  2  or  3,  with  a  maximum  modulation  of  stresses 
of  about  1  MPa.  Release  waves  originating  from  the  lateral  edge  of  the  snow  sample  do  not 
affect  the  stresses  at  the  gauge  plane  positions.  The  propagation  time  to  the  center  axis  for 
release  waves  generated  at  the  snow’s  lateral  edge  is  far  longer  than  the  duration  of  the 
experiment. 
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APPENDIX  B:  Modified  Soil  and  Crushable  Foams  Model 


Preliminary  efforts  to  analyze  our  experimental  measurements  using  the  standard  Soil  and 
Crushable  Foams  model  indicated  that  its  unloading  and  reloading  paths  were  too  restrictive 
(Johnson  et  al.  1990).  Geological  materials  exhibit  nonlinearities  when  they  are  unloaded 
after  shock  wave  passage  (Johnson  and  Green  1976,  Larson  and  Anderson  1979,  Larson 
1 984) .  Additionally,  porous  materials  can  have  a  hysteresis  when  unloaded  and  then  reloaded 
(Johnson  and  Green  1976).  We  modified  the  Soil  and  Crushable  Foams  model  to  include  one 
or  more  linear  or  nonlinear  unloading  curves.  Reloading  can  occur  back  along  the  unloading 
curve  (no  hysteresis)  or  along  chords  that  connect  the  unloading  curve  to  the  /p(£  y3*)  e™  3,1 
state.  On  loading,  the  P-V  curve  is  defined  by  /  p(£yax),  as  was  the  case  for  the  original  model . 

Pressure  is  a  function  of  volumetric  strain  and  is  given  by 


/P(e™ax) 


/p(£T) 


c  _  -max 
tvl  “ 


(  /uj  (€vt) 


eVt<£ir 

EvUi_  ,  <C^VU 
£V(<  £vt_  i 


i  =  1,2,  3, 


f*i  (£v,) 


evt<£r 

1_l<^vaX^vUi 
£vt>  £Vt_, 


where  £v,  is  the  volumetric  strain  at  the  current  time  step  and  eV(  _  j  is  the  volumetric  strain 
at  the  previous  time  step.  The  volumetric  strain  limits  over  which  unloading  fu  t  (e  v,)  and 
reloading  fx ;  (ev,)  functions  are  valid  are  given  by  EVu .  and  e  Vu .  ]  (Fig.  Bl). 

For  loading,  the  P-V  curve  is  given  by 


/tKaX)  =  'M«?aX-evn_i)+  X  S*(evk-evk-l)  ”  =  1.2,3,...  (B2) 


fk-^k-i  *2:1 

iEvk“  £vk-i 


where  and  eVk  are  specified  pressure  and  volumetric  strain  values  defining  the  P-V 
loading  curve  (Fig.  B 1)  and  the  index  n  is  determined  by  the  condition  that 


ev„  - 1  <  £v,  ^  £vn  • 
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Figure  Bl .  Pressure  vs  volumetric  strain  in  terms  of  user-defined  curves 
for  pressure  at  the  maximum  volumetric  strain  [fp^E^1*)] ,  pressure  upon 
unloading  [fU|(ev,)]  and  pressure  upon  reloading[ fr;  (ev<)]  for  tfie 
modified  Soils  and  Crushable  Foams  material  model. 


Tensile  failure  strain  is  initially  set  to  fV/ifO,  where  P(r  is  the  tensile  failure  pressure  in 
the  snow  and  Kq  is  the  snow’s  initial  bulk  modulus.  Failure  strain  is  recalculated  for  each 
increase  in  e(Jiax  by  unloading  from  /p^4*)  using  /„  j(eV|). 

Unloading  P-V  curves  are  defined  by 


i- 1 

/ui  (®vt)  =  j  (e vt "  j- l  )+  ^  ^im(^im) 

m  =  0 


i=  1,2,3,  ... 
j=  1,2,3,  ... 


where 

p 

eij  -  I  -  efr_  S  AEim 

^il  m  =  0 

is  the  strain  at  the  beginning  of  the  curve  segment  containing  eV(  and 


(B4) 


(B5) 


AEj 


0  m  =0 


^im~  F;  m  _  i 


m>0 


(B6) 
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is  the  incremental  strain  across  each  curve  segment  in  the  unloading  curve.  The  specified 
pressure  Pim  and  slope  Kim  define  the  P-V  unloading  curve,  where  Kim  =  0  when  m  =  0. 

The  index  i  specifies  the  unloading  curve  of  interest  and  is  determined  by  the  condition 
that 


£vu  i  _  !  <  E'vaX  -  £v„i 


(B7) 


The  index  j  specifies  the  curve  segment  over  which  a  given  modulus  is  defined  (Fig.  Bl), 
where  j  is  determined  by  the  condition  that 


£j  j_  1  <  £yj  ^  £i  j 


(B8) 


Reloading  can  be  in  the  reverse  direction  along  the  unloading  curve,  for  which  eq  B3  applies 
and  there  is  no  hysteresis.  Alternatively,  reloading  can  be  along  chords  connecting  the 
unloading  curve  to  the  loading  curve,  producing  a  hysteretic  unloading-reloading  cycle.  The 
slope  of  the  reloading  chord  is  required  to  be  greater  than  or  equal  to  the  maximum  slope  of 
/p(£IJJax)between£Vu .  ]  <  e^)4*  <  £Vu  ..Our  analysis  used  reloading  along  chords  where 
the  reloading  curve  is  given  by 

/r  j(® vj  =/r  j  Vt  _  j  )  "*"  S  A£V( 

s  _/f(^aX)-/r,(£v,_  [)  (B9) 

elTiaX  e 
fcV  “  bvt  _  J 


A£yt  —  £v,  —  £v,  _  i  • 

Figure  5  illustrates  the  loading,  unloading  and  reloading  curves  used  in  this  study,  and  Table 
4  gives  the  values  for  P e^,  P,m  and  K;m. 
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